import numpy as np
from scipy import stats
import pandas as pd
'''
公式见书《应用多元统计分析》78、79
'''
dfx = pd.read_csv(r'm_grade.csv')
dfy = pd.read_csv(r'w_grade.csv')

# 转置矩阵x
x = np.mat(dfx.mean()).T
# 转置矩阵y
y = np.mat(dfy.mean()).T

# 男生成绩数据行数
n = dfx.shape[0]
# 女生成绩数据行数
m = dfy.shape[0]
# 协方差tempA1
tempA1 = dfx.cov() * 10
# 协方差tempA2
tempA2 = dfy.cov() * 10

# 中间变量D2
D2 = (n + m - 2) * (x - y).T * np.linalg.inv(tempA1 + tempA2) * (x - y)

# 统计量T2
T2 = (n * m / (n + m)) * D2
# 校验统计量F
F = ((n + m - 4 - 1) / ((n + m - 2) * 4)) * T2
# 获取矩阵中的数据
f = np.array(F)[0]
# 置信水平
a = 0.01
# 临界值
T = stats.f.ppf(1 - a, 4, 17)

# 取小数点后6位
print("统计量=", round(f[0], 6), "临界值=", round(T, 6))
if F > T:
    print("在α=0.01的置信水平下拒绝原假设")
else:
    print("在α=0.01的置信水平下接收原假设")
